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Water  management  is  critical  to  achieving/maintaining  high  performance  of  polymer  electrolyte  fuel  cells 
(PEFCs);  and  elucidating  the  dynamic  behavior  of  liquid  water  droplets  in  a  PEFC  channel  is  essential 
to  water  management.  In  this  work,  the  dynamics  of  liquid  water  droplets  in  a  single  PEFC  gas  flow 
channel  is  investigated  through  theoretical  and  numerical  analyses.  Forces  on  water  droplet,  droplet 
deformation  and  detachment  are  examined.  The  pressure  and  viscous  drags  are  computed  and  compared 
at  different  flow  regimes  (which  exhibit  different  droplet-dynamic  scenarios)  such  as  that  in  the  entrance 
and  fully  developed  flow  regions.  The  expression  for  describing  droplet  shape  change  is  derived,  and  it  is 
found  that  the  droplet  can  deform  significantly  at  high  gas-flow  rates  and  when  the  droplet  is  relatively 
large  (relative  to  the  channel  dimension).  The  detachment  velocity  is  analyzed  by  comparing  the  wall 
adhesion  and  drag  forces,  and  an  expression  relating  the  Weber  number  to  the  Reynolds  number  using 
the  detachment  velocity  is  developed.  Follow-on  work  is  also  briefly  discussed. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Water  management  is  one  of  the  key  issues  in  polymer  elec¬ 
trolyte  fuel  cells  (PEFCs):  appropriate  humidification  is  critical  to 
achieving  high  ionic  conductivity  of  membrane  whereas  excessive 
water  causes  flooding  and  consequently  reduces  cell  performance. 
In  the  gas  flow  channel,  droplets  can  form  at  the  hydrophobic  GDL 
(gas  diffusion  layer)  surfaces  and  hinder  the  transport  of  oxygen 
and  hydrogen  towards  the  respective  catalyst  layers  where  the 
electrochemical  reactions  occur  [1-4]. 

Liquid  water  can  be  produced  in  the  cathode  catalyst  layer  by  the 
oxygen-reduction  electrochemical  reaction  (l02  +  2IT+ +  2e- 
H20)  or  in  the  components  where  condensation  occurs  [5-7]  when 
the  local  temperature  reaches  the  dew  point.  Liquid  water  can 
transport  through  the  porous  media,  which  is  driven  by  the  cap¬ 
illary  force  or  through  the  membrane  by  the  hydraulic  permeation. 
Liquid  water  will  eventually  reach  the  interfaces  between  the 
gas  flow  channels  (GFCs)  and  gas  diffusion  layers  (GDLs),  where 
it  either  forms  droplets  at  the  GDL  surface  or  attaches  to  the 
hydrophilic  channel  wall  [8].  Liquid  water  in  the  GFC  can  block 
local  channels  and  thus  reduce  fuel  cell  performance.  Wang  et  al. 
[9]  indicated  that  the  voltage  fluctuation  occurred  when  the  flow 
stoichiometry  ratio  is  reduced  to  around  2  at  which  the  air  flow 
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rate  is  insufficient  to  remove  liquid  water.  Due  to  the  GDL  sur¬ 
face  roughness  and  contact  angle  hysteresis,  droplets  can  be  held 
by  surface-tension  forces  on  the  GDL  surface,  and  grow  their  bulk 
volume  until  they  are  removed  by  the  gas-phase  reactant  flow 
[10-12].  Several  factors  such  as  the  gas-phase  velocity,  droplet 
diameter,  contact  angle,  contact-angle  hysteresis  and  gas  flow 
properties  can  affect  droplet  detachment.  Droplets  can  be  deformed 
by  the  forces  imposed  by  gas  flows,  which  can  modify  its  pro¬ 
jected  area  and  hence  the  drag  force  [13-16].  The  presence  of 
droplets  increases  the  gas  flow  pressure  drop.  Adroher  and  Wang 
[17]  reported  that  the  two-phase  pressure  amplifier  can  be  as  high 
as  10  at  low  gas  flow  velocity  (<0.5 ms-1).  They  also  presented 
model  analysis  and  compared  their  predictions  with  that  reported 
by  others. 

Several  numerical  and  analytical  approaches  have  been  pro¬ 
posed  to  elucidate  droplet  removal  on  the  GDL  surface.  Chen 
et  al.  [13,14]  provided  simplified  models  for  prediction  of  droplet 
detachment  velocity,  respectively,  in  the  low-Reynolds  number 
(Re)  regime  (in  which  inertia  is  negligible)  and  in  the  inertia- 
dominating  regime  (Re  ~  1000).  The  models  by  Chen  et  al.  account 
for  the  effects  of  static  contact  angle,  contact  angle  hysteresis  and 
gas-phase  velocity  on  detachment  velocity.  Avoidance  of  channel 
lodging  was  also  studied  for  both  cylindrical  and  spherical  droplets. 
Kumbur  et  al.  [  1 5  ]  presented  a  similar  analytical  method  for  droplet 
instability  with  focus  on  the  contact-angle  hysteresis.  The  effects  of 
PTFE  contents  in  MPL  (micro  porous  layer)  and  Reynolds  number 
on  the  contact-angle  hysteresis  were  investigated.  Theodorakakos 
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Nomenclature 

A  area(m2) 

a  coefficient  for  correlation 

b  coefficient  for  correlation 

D,  d  droplet  diameter  (m,  mm) 

Dh  hydraulic  diameter  (m) 

Eo  Eotvos  number 

F  force  (N) 

g  gravitational  acceleration  (ms-2) 

H  channel  height  (m) 

ftd  droplet  height  (m) 

i  direction  vector 

I  unit  tensor 

L  channel  length  (m) 

Le  entrance  length  (m) 

h  normal  vector 

p  pressure  (pa) 

R  local  curvature  (m) 

r  droplet  radius  (m) 

RH  relative  humidity 

U,  u  velocity  (m  s-1 ) 

W  channel  width  (m) 

Wer  Weber  number 

Greek  letters 

a  Volume  fraction 

y  Surface  tension  (N  m-1 ) 

£  area  ratio 

0  contact  angle  (rad,0) 

fi  viscosity  (kg (ms)-1) 

p  density  (kg  m-3) 

Subscript 

a  advancing 

d  droplet 

c  critical 

H  hysteresis 

in  inlet 

press  pressure 

r  receding 

s  static 

visco  viscosity 


et  al.  [16]  studied  water-droplet  deformation  under  cross-flowing 
air  both  experimentally  and  numerically.  In  their  work,  contact 
angles  measured  experimentally  were  compared  to  that  predicted 
by  the  VOF  (Volume  of  Fluid)  method.  He  et  al.  [18]  suggested 
a  simple  analytical  model  to  predict  the  size  of  single  droplet 
for  detachment  under  certain  air  flow  velocity.  Between  advanc¬ 
ing  and  receding  contact  angles,  continuous  contact-angle  change 
along  with  circular  contact  line  was  considered.  The  effect  of  GDL 
hydrophobicity  and  surface  tension  on  detachment  was  studied  by 
He  et  al.  Schillberg  and  Kandlikar  [19]  provided  a  comprehensive 
review  on  analytical  models  for  droplet  detachment  mechanism 
and  suggested  considering  more  accurate  geometry  (droplet  pro¬ 
jected  area,  droplet  surface  area)  for  force  estimation  to  enhance 
the  analytic  models. 

Droplet  formation  occurs  in  the  two-phase  region,  where  the 
gas  flow  reaches  the  saturated  state.  Wang  [8]  gave  the  prediction 
of  the  onset  of  two-phase  flow  when  dry  air  is  used.  Droplets  tend 
to  form  in  preferential  sites  after  the  flowing  gas  mixture  reaching 
the  saturated  state  and  they  are  not  uniformly  distributed  along 
the  channel  [20,21].  Ous  and  Arcoumamis  [21]  studied  the  effect 


of  air  velocity,  relative  humidity  and  external  load  of  fuel  cells  in 
the  presence  and  removal  of  water  droplet  through  visualization. 
Droplets  are  formed  non-uniformly  along  the  channel  and  air  veloc¬ 
ity  is  inversely  proportional  to  the  detachment  droplet  dimension. 
Shirani  and  Masoomi  [22]  studied  the  main  factors  of  2-D  droplet 
deformation  by  using  the  VOF-PCIL  method.  It  is  found  that  droplet 
deformation  in  the  air  flow  direction  is  highly  related  to  the  capil¬ 
lary  number.  Zhu  et  al.  [23,24]  used  the  2-D  and  3-D  VOF  methods 
to  investigate  the  impact  of  pore  and  channel  geometry  on  droplet 
dynamics.  Fang  et  al.  [25]  emphasized  the  importance  of  contact- 
angle  hysteresis  on  detachment  velocity,  showing  that  the  VOF 
without  taking  into  account  of  contact-angle  hysteresis  can  cause 
low  detachment  velocity. 

Despite  previous  efforts,  a  comprehensive  study  on  droplet 
dynamics  that  delineates  the  forces  over  a  droplet  at  different 
regions  of  channel  flows,  and  theoretical  analysis  on  droplet  defor¬ 
mation  and  detachment  is  highly  needed.  In  this  work,  we  examine 
the  forces  over  a  droplet  at  various  conditions  and  different  loca¬ 
tions.  Analysis  is  performed  to  predict  water  droplet  deformation, 
and  numerical  models  of  droplet  detachment  are  also  introduced. 
The  influence  of  droplet  diameter,  gas-phase  velocity,  and  location 
of  liquid  droplet  are  presented  for  both  air  and  hydrogen  flow  fields. 
We  will  compare  the  analytical  solutions  of  the  droplet  deforma¬ 
tion  and  removal  with  experimental  data  and  the  VOF  simulation 
prediction  in  a  sequel  paper. 


2.  Droplet  dynamics 

In  this  section,  three  major  topics  of  droplet  dynamics  are 
described  and  discussed,  which  are  important  to  elucidating 
droplet  behaviors  in  PEFC  gas  flow  channels. 


2A.  Forces  on  a  liquid  water  droplet 


A  water  droplet  on  the  GDL  surface  in  the  gas  flow  channel 
experiences  drag  forces  exerted  by  the  surrounding  gas  flow.  Since 
channel  flow  is  viscous,  the  sharp  change  of  the  gas  velocity  near 
the  droplet  surface  causes  a  drag  force  over  the  droplet.  Drag  force 
can  also  arise  from  the  gas  pressure:  the  upstream  has  a  higher 
pressure,  leading  to  the  pressure  drag.  In  the  gas-flow  channel  of  a 
PEFC,  the  surface  tension  is  usually  a  dominant  factor  that  controls 
the  liquid-droplet  shape.  Adamson  and  Gast  [27]  examined  the  free 
droplet  shape  considering  surface  tension  and  gravitational  force 
in  static  air.  The  effect  of  the  gravitational  force  is  usually  small  for 
droplets  with  diameters  on  the  order  of  a  few  hundred  microns  in  a 
millimeter-sized  channel,  which  can  be  explained  using  the  Eotvos 
number  defined  as: 


Y 


(1) 


The  cases  considered  in  the  present  work  show  a  relatively 
low  Eotvos  number  (Eo:  0.001 5-0.1 5)  with  which  the  gravitational 
effect  is  negligible  [30].  Consequently,  the  droplet  shape  in  slow  air 
flow  can  be  assumed,  to  a  good  approximation,  to  be  spherical  and 
has  a  constant  curvature  at  any  local  point.  Fig.  1  shows  the  force 
balance  on  a  spherical  liquid  droplet.  At  the  onset  of  detachment 
(that  is,  immediately  before  the  droplet  is  detached),  all  forces  on 
a  droplet,  including  the  wall  adhesion,  are  balanced  to  keep  the 
droplet  stationary.  The  drag,  created  by  air  flows  near  the  droplet, 
is  the  sum  of  the  pressure  and  viscous  forces.  The  pressure  force 
can  be  calculated  by: 


Fpress  — 


p  dA 


(2) 
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Fig.  1.  Force  balance  on  a  liquid  water  droplet  on  a  hydrophobic  wall. 


where  i  is  the  direction  vector  of  gas  flow,  n  the  normal  vector  of 
the  infinitesimal  surface,  p  local  pressure  and  A  the  area  of  droplet 
surface.  By  employing  CFD  (computational  fluid  dynamics)  with  a 
fine  mesh,  the  pressure  distribution  on  the  droplet  surface  can  be 
integrated  over  the  droplet  surface  in  the  flow  direction  to  com¬ 
pute  the  drag  force  caused  by  pressure.  Another  force,  the  viscous 
friction  at  the  gas-liquid  interface,  can  be  evaluated  by: 


l 


Mgas  V u  dA 


(3) 


The  local  velocity  gradient  is  calculated  based  on  the  CFD  results 
of  the  velocity  field.  The  total  drag  is  then  the  summation  of  pres¬ 
sure  force  and  viscous  force: 


^drag  —  Fpress  +  ^Visco  (4) 

Wall  adhesion  (or  force  arising  from  surface  tension)  can  be 
divided  into  the  normal  (z-direction)  and  tangential  (x-direction) 
components.  The  maximum  tangential  force  is  determined  by  the 
static  contact  angle  (0S),  contact-angle  hysteresis  (0H)  and  surface 
tension  (y)  [13]: 

.  ?  .  0U 

Fyx  —  Fyx, r  +  Fyx, a  =  2yjtd  •  sin  0S  •  sin  —  (5) 

For  the  hydrophobic  GDL  surface,  the  effects  of  advancing  and 
receding  angles  are  opposite  to  each  other  and  thus  partly  canceled 
out.  When  there  is  surrounding  gas  flow,  the  advancing  angle  is 
greater  than  receding  angle,  making  the  combined  wall  adhesion 
opposite  to  the  gas  flow  direction.  The  normal  wall  adhesion  can  be 
estimated  by  integrating  the  adhesion  force  over  the  contact  line 
of  the  droplet  and  wall: 


Fyz  =  2y-  sin  0S^  •  J  sin  ^0r  +  (p^j  d(f  (6) 

0  is  the  angle  of  the  projected  circle  on  the  GDL  surface.  In  the 
present  work,  a  simple  relation  is  adopted  between  the  static  con¬ 
tact  angle  and  contact-angle  hysteresis: 

er  =  Os  -  (7) 

A  local  contact  angle  f  can  be  defined  along  the  contact  line  on 
the  GDL  surface  as  [15]: 


$  =  er  + 


9,-9^ 


<t> 


(8) 


Eq.  (6)  can  be  re-organized  as: 

1 

71 

=  2 yf-d  ■  sin  6S  ■  sin  ^  •  sin  ^  (9) 

C7h  2  2 


This  force  component  is  the  maximum  that  can  balance  the  lift 
force  over  the  droplet  due  to  the  surrounding  flow.  When  the  lift 
force  overcomes  this  barrier,  the  droplet  will  be  detached  from  the 
surface.  The  criteria  for  water  droplet  detachment  can  be  defined 
by  considering  two  directions: 


^drag  >  Fyx 


(10) 


^Lift  >  Fyz 


(11) 


Satisfaction  of  either  of  the  two  criteria  will  lead  to  droplet 
detachment.  In  reality,  the  second  criterion  is  more  difficult  to 
satisfy  than  the  first  one.  Thus,  the  first  criterion  is  frequently 
used  for  droplet  detachment  analysis  and  the  critical  condition,  i.e. 
Fdrag  =  Fyx,  can  be  expressed  by: 


?■  f(n 

JA 


?  0H 

P  +  Mgas  Vu)dA  =  2ynd  •  sin  0S  •  sin  — 


(12) 


2.2.  Droplet  deformation  due  to  pressure  variation 


Droplets  can  be  deformed  by  the  forces  exerted  by  the  surround¬ 
ing  air  flow,  causing  the  droplet  to  deviate  away  from  the  spherical 
shape  [22-25].  There  are  two  major  forces  by  air  flow:  one  is  the 
viscous  force,  the  other  is  pressure  variation.  The  former  is  due 
to  the  viscosity  of  air  (hydrated  with  liquid  water),  which  causes 
the  sharp  change  of  velocity  near  the  droplet  surface,  and  is  usu¬ 
ally  along  the  tangential  direction  at  the  surface.  This  force  highly 
depends  on  the  local  flow  field  and  can  vary  greatly  when  air  flow 
condition  changes.  At  a  high  air  flow  rate,  backward  air  flows  can 
develop  near  the  liquid  surface.  The  other  force  results  from  the 
droplet  obstruction  which  causes  the  pressure  variation  over  the 
droplet  surface.  The  pressure  difference  at  two  typical  sites  can  be 
estimated:  one  is  at  the  droplet  front  where  the  flow  is  stagnated; 
the  other  is  at  the  side  where  the  air  flow  cross  section  reaches  a 
local  minimum.  The  air  pressure  will  also  change  in  correspond¬ 
ing  to  the  velocity  change.  The  pressure  change  will  alter  the  local 
curvature  of  the  droplet.  Fig.  2  shows  schematic  of  a  liquid  droplet 
in  a  gas  flow  channel.  The  local  curvature  of  the  droplet  can  be 
calculated  by  the  Young-Laplace  Equation: 

Pliq  -  Pgas  =  (13) 

Assuming  a  constant  liquid  pressure  throughout  the  droplet  at 
steady  state,  Eq.  (13)  can  be  differentiated  to  evaluate  the  impact 
of  the  gas  pressure  variation  dpgas  on  the  local  curvature: 


dR 


R 

2y 


dpgas 


(14) 


where  R  is  the  local  curvature  of  the  gas-liquid  interface.  Evaluation 
of  dpgas  can  follow  that  in  the  orifice  plates  or  Venturi  tubes  [26]: 
based  on  the  Bernoulli’s  principle,  the  gas-phase  pressure  drop  is 
proportional  to  the  square  of  air  velocity  and  the  local  averaged  gas 
velocity  can  be  computed  using  the  mass  conservation.  To  better 
elucidate  the  issue,  we  also  consider  a  2-D  case,  which  assumes  the 
channel  only  has  finite  dimensions  in  the  along-channel  and  height 
directions,  and  a  2-D  droplet  has  a  shape  similar  to  a  cylinder.  The 
area  ratio  (f )  between  upstream  and  on  the  droplet  for  2-D  and  3-D 


71 


122 


S.C.  Cho  et  al.  /  Journal  of  Power  Sources  206  (2012)  1 19-128 


Fig.  2.  Schematic  of  a  liquid  water  droplet  in  a  2-D  and  3-D  micro  gas  channel. 


can  be  evaluated  by: 

a, 


2-D  :  fto  =  t1  = 


WxH 


A2  W  x[H  -  R0(l  —  cos  0S)] 

H 


H  -  R0(l  —  cos  9S 


3-D  :  f3D  =  v1  = 


WxH 


A3  WxH  —  R2{0s  -  cos  Os  sin  0S) 


(15) 


(16) 


where  Ro  is  the  radius  of  curvature  without  deformation.  Substitu¬ 
tion  of  this  area  ratio  into  the  Bernoulli’s  equation  gives: 


2- D:  dpgas  =  P1-P2  =  ^22D 

3- D:  dpgas=P1-P3  =  ^WfD 


Substituting  the  pressure  differences  into  Eq.  (14)  yields: 

2D-  —  =  —t2  =  —  t2 

‘  R  4 y  V2D  4  fV2D 

dR  pujRo  R  2  Wer  R  2 
'  R  4y  R0?3D  4  iV3D 


(17) 

(18) 

(19) 

(20) 


where  Wer  is  the  Weber  number,  defined  as  the  ratio  between  fluid 
inertia  and  surface  tension: 


Wer  = 


pu2R o 


(21) 


2.3.  Droplet  detachment 


The  gas  pressure  drop  across  droplet  can  be  related  to  the  over¬ 
all  pressure  drop  by  using  average  velocity,  viscosity  and  channel 
height  as  [13] 


(p'o  -  P'L)  =  (Po  -  Pl)  -  (L  -  l)^-  (24) 


(P'o  -P'0=^>  (25) 

Based  on  the  mass  balance  (LT  =  (B/b)L0,  Eq.  (24)  can  be 
expressed  as  [13] 

(Po-P'l)  =  (Po-Pl) - - - 3  (26) 

l  +  (L  —  l)(b/B)3 

Substitution  of  Eq.  (26)  into  (23)  yields  [13] 
r  ,  ,  l2(2  B-b) 

-Fdrag-(Po-PL)(  +  (L_()(b/B)3  (27) 

By  defining  dimensionless  droplet  height  as  H  =  r(1- 
cos  0a)/2B,  Eq.  (27)  becomes: 


(28) 


_F  (Po-PlWA)  I(l+H) 

drag  (1  —  H)3  +  (//L)[l  —  (1  —  H)3] 

For  a  sufficiently  large  droplet  with  0a  >90°,  l  =  2r,  the  viscous 
drag  can  be  further  simplified  in  terms  of  average  velocity  in  the 
upstream  and  downstream  regions  as  [13] 


48  fiUB 


(1  +  H)H2 


p,  — 

ras  (1  -cos6>a)(i  _C0SPa)(l  -H)3+(4B/L)H[1  -(1  -H)3] 

(29) 

The  adhesion  force  on  the  contact  line  between  the  wall  and 
droplet  can  be  expressed  by  considering  the  advancing  and  reced¬ 
ing  contacts: 


2.3.  t.  Analysis  using  the  control-volume  approach 

In  the  occasions  that  the  droplet  shape  change  is  negligible, 
the  detachment  velocity,  defined  as  the  minimum  gas  velocity  to 
detach  the  droplet,  i.e.  Fdrag  =  F)/x,  can  be  estimated  based  on  the 
force  balance  over  a  spherical  droplet.  Following  Chen  et  al.  [13], 
the  force  balance  along  the  channel-flow  direction  of  the  control 
volume  can  be  expressed  as  (see  Fig.  3) 

(Po  -  P'l)2B(  +  l2T™  +  Fdrag  =  0  (22) 

where  l2  x™  and  Fdrag  represent  the  shear  stress  on  the  top  wall  and 
drag  on  the  water  droplet.  By  assuming  fully  developed  laminar 
Newtonian  gas  flow,  the  viscous  drag  can  be  approximated  as  [13] 

(23) 


Fyx  =  7r(r  sin  6a)y  cos(180  -  0a)  +  7t{r  sin  6a)y  cos  0r  (30) 

When  the  surface  tension  force  is  larger  than  the  viscous  drag 
balance,  the  droplet  will  be  detached,  therefore  the  critical  velocity 
can  be  calculated  by  setting: 


Fdrag  —  F^  (31) 

By  defining  y  =  sin(l/2)(0a  -0T),  substituting  Eqs.  (29)  and  (30) 
into  (31)  yields: 


y \/ 1  - y2  -  (cot  6s)y2  — 

Tty  sin  9S 

_ tf(l+H) _ 

(1  -  cos  0S)(1  -ff)3  +(4B/I)H[1  -(1  -ff)3] 


(32) 


Fdrag  —  (Po  Pl)(2®  ^)( 
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Fig.  3.  Control  volume  enclosing  the  liquid  water  droplet  [13]. 


This  above  expression  can  be  further  rearranged: 


3.  Numerical  simulation 


yV i  -y2  -  (cot  0s)y2  -  Re  1  Wer — — — 

r  7T  sin  0S 

_ H(l+H) _ 

(1  -  cos  0S)(1  -  H)3  +  (4£/L)H[l  -  (1  -  ff)3] 


(33) 


2.3.2.  Derivation  using  the  drag  coefficient  (CD ) 

Another  method  that  can  be  used  to  estimate  the  droplet  detach¬ 
ment  velocity  the  drag  coefficient  CD  defined  by: 


_  ^drag/^P 

“  (l/2)pl/2 


(34) 


where  Ap  is  the  projected  area  of  the  droplet  which  is  estimated 
as  Ap  =  (0S  -  sin  0S  cos  0s)d2l 4.  The  drag  coefficient  is  determined  by 
the  Reynolds  number.  Chen  [14]  adopted  the  following  formula: 


Q> 


30 

V^h 


ReH  = 


pUH 

li 


(35) 


The  surface  tension  force  (Eq.  (30))  can  be  simplified  as: 


Fyx  =  Tcdy  sin3  0S  sin  -(0a  -  0r 


(36) 


In  order  to  evaluate  the  forces  on  a  liquid  droplet  in  a  gas  flow 
channel,  the  3-D  Navier-Stokes  equations  for  incompressible  vis¬ 
cous  flow  are  solved  numerically: 


V  •  (pv)  =  0 


(39) 


V  •  (pvv)  =  -Vp  +  V  •  /x 


(Vu  +  Vur)  -  -V  •  vl 


+  pg 


(40) 


In  this  evaluation,  we  assume  negligible  deformation,  thus  the 
droplet  has  a  spherical-cap  shape.  We  also  assume  the  non-slip 
condition  of  the  gas  flow  on  the  droplet  surface.  Fig.  4  shows 
the  computational  domain  with  boundaries.  The  computational 
domains  are  built  with  about  1.3  million  hexahedral  cells  (which 
was  arrived  at  after  a  grid-independent  study).  Air  flow  is  intro¬ 
duced  through  uniform  velocity  at  inlet  boundary  uniformly  and 
it  exits  through  pressure  outlet  boundary.  Due  to  symmetry,  only 
half  of  the  channel  is  simulated.  Non-slip  condition  is  applied  to 
all  wall  boundaries.  Various  droplet  sizes  are  considered  and  air 
flow  velocity  is  set  between  0.5  and  10.0  ms-1  based  on  previous 
studies  (e.g.  [13,14]).  Details  of  the  cases  considered  for  numerical 
simulation  are  presented  in  Table  1 .  A  commercial  CFD  code  (Fluent 
6.3.26)  was  used  to  solve  the  highly  nonlinear  governing  equations. 


4.  Results  and  discussion 


Again,  the  critical  gas  velocity  can  be  obtained  by  setting 

^drag  =  Fyx  [14]: 


Uc=  — 


'  H  - 

1/3 

4  ny  sin2  0S  sin(l/2)(0a  -0r) 

-PI1  - 

15  (0S  -  sin  0S  cos  0s)d 

2/3 


(37) 


A  simplified  relation  between  the  Weber  number  and  Reynolds 
number  can  be  obtained  from  the  above  equation: 


w P  dp-0.5  _  2  jrsin  e_s  sin(l/2)(6>a-6>r) 
r  H  15  (9S  -  sin  9S  cos  9S) 


(38) 


4. 1 .  Force  components 

In  real-world  PEFC  operations,  the  maximum  possible  air  veloc¬ 
ity  in  a  channel  with  a  static  droplet  on  the  GDL  surface  is  limited 
by  the  detachment  velocity.  For  a  parametric  study  on  the  forces 
over  a  droplet,  the  gas  velocity  is  allowed  to  exceed  that  limit 
while  the  droplet  remains  static.  Fig.  5  shows  the  two  drag  forces 
(Eqs.  (2)-(4))  as  a  function  of  the  gas  velocity.  Various  droplet  sizes 
are  compared  and  results  are  presented  for  fully  developed  region 
(Case  1)  and  the  entrance  region  (Case  2).  As  expected,  both  pres¬ 
sure  and  viscous  force  increase  rapidly  with  the  gas  velocity.  For 
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Table  1 

List  of  numerical  simulation  cases. 


Case  no. 

Working  fluid 

Flow  condition 

Channel  dimension 

(mm3) 

Contact  angle 

Contact  angle 
hysteresis 

Droplet  diameter 
(mm) 

Air  flow  velocity 
(ms-1) 

1 

Air 

Fully  developed 

2 

3 

Air 

Hydrogen 

Entrance  length 
Fully  developed 

1.8  x  1.0x23 

150° 

10° 

0.1-0.6 

0.5-10.0 

0.0000008 

0.0000007  - 

0.0000006  - 

0.0000005  - 

^  0.0000004  - 

8 

£  0.0000003- 
0.0000002  - 

0.0000001  - 
0.0000000  - 

- - 1 - - - 1 - - - r 

0  2  4  6 

Air  velocity  [m/s] 


0.0000030  - 


0.0000025  - 


0.0000020  - 


0.0000015- 


0.0000010 


0.0000005  - 


0.0000000  - 


Entrance  Length  -  Total 
Fully  Developed  -  Total 
Entrance  Length  -  Viscous 
-O—  Fully  Developed  -  Viscous 

☆  Entrance  Length  -  Pressure 

☆  Fully  Developed  -  Pressure 

D=0.4  mm 


0.000008  -| 
0.000007  - 
0.000006  - 
0.000005  - 
0.000004  - 

8 

£  0.000003- 
0.000002  - 

0.000001  - 
0.000000  - 


- Entrance  Length  -  Total 

- Fully  Developed  -  Total 

-O-  Entrance  Length  -  Viscous 
-O-  Fully  Developed  -  Viscous 

☆  Entrance  Length  -  Pressure 

☆  Fully  Developed  -  Pressure 


D=0.6  mm 


Air  velocity  [m/s] 


Air  velocity  [m/s] 


Fig.  5.  Comparison  of  drag  forces  between  the  entrance  and  fully  developed  regions  (Case  1  and  Case  2). 


S.C.  Cho  et  al. /Journal  of  Power  Sources  206  (2012)  119-128 


125 


(a)  Entrance  Length  (d=0.6mm) 


(b)  Entrance  Length  (d=0.1mm) 


(c)  Fully  Developed  (d=0.6mm)  (d)  Fully  Developed  (d=0.1mm) 


Fig.  6.  Velocity  contours  in  the  entrance  (a  and  b)  and  fully  developed  (c  and  d)  regions,  respectively  (uin  =  3.0  m  s-1 ). 


small  droplets,  the  viscous  force  dominates,  while  the  pressure 
drag  becomes  more  important  for  relatively  large  droplets.  Fig.  5 
also  compares  the  forces  for  two  typical  locations  in  a  channel,  i.e. 
the  entrance  (Case  2)  and  fully  developed  regions  (Case  1),  which 


show  flow  in  distinct  regimes.  Many  studies  only  considered  a  fully 
developed  channel  flow  in  fuel  cells  because  of  the  large  ratio  of  the 
channel  length  to  cross-section  dimension.  In  the  entrance  length, 
the  boundary  layer  is  undergoing  development,  and  due  to  the 


Gaseous  phase  velocity  (m/s) 


Fig.  7.  Drag  forces  for  air  and  hydrogen  flows,  respectively,  at  various  droplet  diameters  (Case  1  and  Case  3). 
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(a)  Effect  of  droplet  height  at  various  gas  velocities  (b)  Effect  of  gas  velocity  for  various  droplet  heights 

Fig.  8.  Droplet  deformation  in  air  flow  (2-D  and  3-D)  calculated  by  Eqs.  (19)  and  (20).  (a)  Effect  of  droplet  height  at  various  gas  velocities;  (b)  effect  of  gas  velocity  for  various 
droplet  heights. 


momentum  diffusion,  the  boundary  layer  will  grow  and  eventu¬ 
ally  meet  at  the  centerline.  After  that,  the  channel  flow  reaches  the 
fully  developed  stage  with  a  parabolic  velocity  profile.  In  PEFCs, 
the  Reynolds  number  of  channel  flows  is  usually  relatively  small 
or  moderate  (that  is,  not  exceedingly  large)  such  that  flow  is  in 
the  laminar  regime.  The  hydrodynamic  entrance  length  IE  can  be 
evaluated  using  that  developed  in  tube  laminar  flow:  LE/d  ~  0.06Re, 
where  d  is  the  tube  diameter.  A  Re  of  200  will  lead  to  LE~12d, 
which  is  small  but  not  negligible.  For  a  1.0  mm  diameter  tube,  the 
entrance  length  is  around  12  mm,  which  is  about  10%  of  a  10  cm 
length  channel.  The  velocity  profiles  near  the  water  droplet  for 
both  near  entrance  and  fully  developed  conditions  are  plotted  in 
Fig.  6.  In  the  entrance  region,  a  small  droplet  (d<0.2mm)  experi¬ 
ences  a  stronger  force  than  that  in  the  fully  developed  region.  This 
can  be  explained  by  that  in  the  entrance  region  the  gas  velocity  is 
almost  uniform  in  the  gas  channel,  therefore  the  velocity  around 
small  droplets  on  the  wall  undergoes  larger  surrounding  gas  flows 
in  comparison  with  the  fully  developed  region.  As  a  result,  the 
detachment  velocity  of  a  small  droplet  is  lower  in  the  entrance 
region  than  that  in  the  fully  developed  condition.  As  the  droplet 
size  increases,  the  pressure  force  becomes  dominant.  In  the  fully 
developed  region,  the  local  air  velocity  in  the  core  of  the  channel 


is  higher  than  the  average  value  that  approximately  equals  to  the 
nearly  uniform  velocity  in  the  entrance  region.  A  large  droplet  will 
physically  extend  to  the  core  region,  facing  a  higher  air  velocity 
in  the  fully  developed  region.  Since  the  pressure  force  is  approxi¬ 
mately  proportional  to  the  square  of  the  gas  velocity  magnitude,  a 
larger  droplet  experiences  a  larger  drag  force,  i.e.  a  lower  detach¬ 
ment  velocity  is  in  the  fully-developed  region  than  in  the  entrance 
region.  When  the  working  fluid  is  switched  to  hydrogen  (i.e.  in  the 
anode  gas  flow  channel),  the  magnitude  of  the  forces  decreases 
by  about  an  order,  as  shown  in  Fig.  7,  due  to  the  lower  gas  den¬ 
sity.  Note  that  many  studies  indicated  that  liquid  can  emerge  on 
the  anode  side  [12,28,29].  If  droplets  eventually  form  on  the  anode 
GDL  surface,  it  can  be  difficult  to  remove  them. 

4.2.  Droplet  shape  change 

To  examine  droplet  shape  changes  under  various  conditions,  the 
ratio  of  local  curvature  change  to  the  initial  surface  curvature  (Eqs. 
(19)  and  (20))  is  plotted  in  Fig.  8.  The  droplet  height  is  normalized 
by  the  channel  height  for  2-D  and  droplet  diameter  is  normalized 
by  the  hydraulic  diameter  of  the  gas  channel  for  3-D.  The  droplet 
deformation  rapidly  increases  with  the  droplet  height,  particularly 
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Fig.  9.  Droplet  deformation  in  hydrogen  channel  flow  (3-D)  calculated  by  Eq.  (20).  (a)  Effect  of  droplet  height  change;  (b)  effect  of  gas  velocity  change. 


when  approaching  the  channel  height.  The  reason  is  that  the  accel¬ 
eration  of  gas  flow  around  the  droplet  due  to  its  obstruction  lowers 
the  local  gas  pressure.  The  increasing  rate  is  much  faster  for  the  2- 
D  case  than  the  3-D  case  because  the  domain  of  the  latter  has  four 
corners  which  can  allow  gas  flow  to  pass  through.  In  addition,  the 
deformation  magnitude  increases  with  the  air  velocity,  as  expected. 
At  low  air  velocity,  the  forces  by  the  gas  flow  are  weak,  leading 
to  a  small  Weber  number  and  hence  small  deformation.  It  can  be 
seen  that  the  deformation  is  usually  smaller  than  1 0%  when  the  gas 
flow  velocity  is  less  than  1 .0  m  s-1  for  both  2-D  and  3-D  cases.  Even 
for  the  cases  with  velocity  smaller  than  5ms-1,  the  droplet  defor¬ 
mation  is  smaller  than  10%  for  the  3-D  cases.  Therefore,  to  a  good 
approximation  we  can  neglect  the  effect  of  droplet  shape  change 
in  gas  flow  channels  for  most  operating  condition,  which  simplifies 
the  analysis.  Fig.  9  shows  the  parametric  study  of  droplet  deforma¬ 
tion  in  hydrogen  gas  flow.  The  magnitude  of  droplet  deformation 
is  smaller  than  that  for  the  flowing  air,  which  is  due  to  the  smaller 
density  of  hydrogen  flow. 


By  using  this  expression,  Eq.  (38)  can  be  changed  to: 


u  4  7T 

WerRef,  = - 

H  a 


sin2  0S  sin(l /2)(0a  -  0r) 
(0S  -  sin  Os  cos  0S) 


(42) 


Fig.  1 1  plots  the  detachment  velocity  for  various  droplet  sizes 
for  both  air  and  hydrogen  flows  using  Eq.  (42).  It  is  observed  that 
the  detachment  velocity  dramatically  increases  as  droplet  diame¬ 
ter  decreases.  This  trend  is  well  matched  to  the  previous  findings 
reported  in  the  literature  [16,17].  The  detachment  velocity  is  also 
sensitive  to  the  contact  angle  and  contact-angle  hysteresis:  as  the 
contact  angle  decreases  and  the  contact-angle  hysteresis  increases, 
the  detachment  velocity  becomes  larger.  Since  the  hydrogen  flow 
exerts  a  smaller  drag  force  on  a  droplet,  higher  gas  velocities  are 
required  to  detach  water  droplet  in  the  anode  side.  Fig.  12  dis¬ 
plays  the  comparison  of  the  Re-We  relations  at  the  detachment 
velocity  at  various  droplet  sizes  among  the  models,  i.e.  Eqs.  (33), 
(38)  and  (42).  Contact  angle  and  contact-angle  hysteresis  are  set  at 
1 50°  and  1 0°,  respectively,  and  channel  length  is  assumed  as  4.5  cm 


4.3.  Detachment  velocity 

The  drag  coefficient  CD  is  a  dimensionless  parameter  character¬ 
izing  the  drag  force  on  an  object  in  fluid  flow.  Many  correlations  of 
the  drag  coefficient  as  a  function  of  the  Reynolds  number  (Re)  have 
been  developed  for  objects  in  unconfined  flows.  Similar  approach 
can  be  taken  for  the  drag  coefficient  over  a  droplet  inside  channels. 
In  the  range  of  Re  for  typical  fuel  cell  operation,  the  drag  coefficient 
of  various  droplet  diameters  are  numerically  calculated  using  the 
3-D  simulation  results,  which  are  plotted  in  Fig.  10.  It  can  be  seen 
that  the  one  used  by  Chen  [14]  (Eq.  (38))  gives  a  fair  agreement, 
but  exhibits  deviation  for  small  droplets  (d~0.1  mm).  It  might  be 
impossible  to  use  one  formula  to  describe  the  full  ranges  of  gas  flow 
and  channel  dimension.  A  generalized  approach  can  be  expressed: 

CD  =  a  x  Re^  (41) 

In  this  study,  the  fitting  parameters  for  the  fully  developed  con¬ 
dition  are: 

°-1757  f  D\ 

,  b  =  0.2158  x  l-j  -0.6384 


a  =  46.247  x  (R) 


Fig.  10.  Drag  coefficients  in  a  micro  channel  at  various  conditions. 
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Fig.  11.  Detachment  velocities  in  air  flow  and  hydrogen  (H2)  flows,  respectively, 
calculated  by  Eq.  (42). 


Fig.  12.  The  We-Re  relation  at  various  droplet  sizes  calculated  by  Eqs.  (33),  (38)  and 
(42),  respectively. 


for  comparison.  It  can  be  seen  that  the  Weber  number  increases 
with  the  Reynolds  number  at  the  detachment  velocities.  As  droplet 
size  increases,  the  Weber  number  decreases  since  less  drag  force  is 
required  for  detachment.  The  correlation  used  by  Chen  [14]  shows  a 
higher  Weber  number,  which  gives  the  upper  bound  of  the  droplet 
stable  region. 

While  several  important  analytical  solutions  of  droplet  dynam¬ 
ics  such  as  deformation  and  detachment  are  obtained  in  the  present 
work,  it  is  important  to  compare  model  prediction  with  experimen¬ 
tal  data.  In  the  follow-on  work  as  presented  in  our  sequel  paper, 
we  will  present  extensive  experimental  work  on  droplet  deforma¬ 
tion  and  detachment  and  comparison  between  experimental  and 
simulation  (or  model  prediction)  results.  In  addition,  the  VOF  sim¬ 
ulation  provides  detail  of  two-phase  flow  along  with  the  interface 
track.  We  will  report  results  from  both  2-D  and  3-D  VOF  studies  and 


compare  their  predictions  with  the  analytical  solutions  presented 
in  this  paper. 

5.  Conclusion 

In  this  study,  water  droplet  dynamics  in  a  PEFC  gas  flow  chan¬ 
nel  was  analyzed.  The  drag  forces  on  a  droplet  were  computed 
and  compared  under  various  conditions.  We  found  that  the  vis¬ 
cous  force  has  significant  impact  on  a  small  droplet  in  low  gas 
velocity  regime  whereas  the  pressure  drag  is  important  for  a  large 
droplet  when  the  gas  flows  relatively  slowly.  Small  droplets  in  the 
entrance  region  get  more  drag  than  that  under  fully  developed  con¬ 
dition  while  a  large  droplet  experiences  greater  drag  under  the  fully 
developed  condition.  When  hydrogen  is  used  as  the  flowing  gas, 
both  pressure  and  viscous  drags  are  reduced.  Droplet  deformation 
was  also  analyzed  and  a  formula  was  derived  to  quantify  the  defor¬ 
mation  caused  by  pressure  variation  at  two  typical  locations.  We 
found  that  droplet  deformation  is  fairly  small  in  most  cases,  but 
significant  shape  change  may  occur  at  high  gas  velocities  and  large 
droplet  sizes.  Further,  the  drag  coefficient  over  a  droplet  in  a  gas 
flow  channel  was  computed  and  a  model  was  developed  to  pre¬ 
dict  the  detachment  velocity.  We  found  that  the  Weber  number 
increases  with  Reynolds  number  at  the  detachment  velocity  and 
decreases  with  growing  droplet  size.  In  the  follow-on  work  that  will 
be  reported  in  a  sequel  paper,  we  will  present  experimental  results 
and  extensive  VOF  simulation  to  validate  the  analytical  solutions 
of  the  droplet  deformation  and  removal  reported  in  the  present 
work. 
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